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Abstract 

The particles model, the collision model, the polynomial development used for the 
equilibrium distribution, the time discretization and the velocity discretization are 
factors that let the lattice Boltzmann framework (LBM) far away from its conceptual 
support: the continuous Boltzmann equation (BE). Most collision models are based 
on the BGK, single parameter, relaxation-term leading to constant Prandtl num- 
bers. The polynomial expansion used for the equilibrium distribution introduces an 
upper-bound in the local macroscopic speed. Most widely used time discretization 
procedures give an explicit numerical scheme with second-order time step errors. 
In thermal problems, quadrature did not succeed in giving discrete velocity sets 
able to generate multi-speed regular lattices. All these problems, greatly, difficult 
the numerical simulation of LBM based algorithms. In present work, the systematic 
derivation of lattice-Boltzmann models from the continuous Boltzmann equation 
is discussed. The collision term in the linearized Boltzmann equation is modeled 
by expanding the distribution function in Hermite tensors. Thermohydrodynamic 
macroscopic equations are correctly retrieved with a second-order model. Velocity 
discretization is the most critical step in establishing regular-lattices framework. In 
the quadrature process, it is shown that the integrating variable has an important 
role in defining the equilibrium distribution and the lattice-Boltzmann model, lead- 
ing, alternatively, to temperature dependent velocities (TDV) and to temperature 
dependent weights (TDW) lattice-Boltzmann models. 
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Introduction 



Following Lallemand & Luo, [1], the presently known lattice-Boltzmann equa- 
tion (LBE) has not been able to handle realistic thermal (and fully compress- 
ible) fluids with satisfaction. Simulation of thermal lattice-Boltzmann equa- 
tion is hampered by numerical instabilities when the local velocity increases. 
Readers are referred to this work for an excellent review of known works on 
thermal and compressible lattice-Boltzmann schemes. 

Rigorously, fluid flow is never isothermal. Consider, for instance, a channel 
Poiseuille flow between two adiabatic solid surfaces. Due to the viscous con- 
version of mechanical in thermal energy, temperature will vary, attaining a 
minimum at the channel symmetry axis, where the local speed is a maximum. 
This temperature variation can be very small, but increases with the local 
average speed u and with the fluid viscosity. If it is agreed that the tempera- 
ture in a given site is related to the expected value E of the lattice-particles 
fluctuation kinetic energy, this temperature variation is drafted in LB ather- 
mal simulation, since E varies from site to site in accordance with the local 
macroscopic speed, u, attaining a minimum where u is maximum. However, 
in athermal lattice-Boltzmann models, this thermal spatial non-equilibrium 
is not properly compensated by heat flow because athermal models were not 
conceived for correctly describing energy transfer. In this manner, since tem- 
perature gradients cannot be avoided in athermal LB simulation, in contrast 
with classical CFD isothermal simulation, they may become sources of numer- 
ical instability. 

In conclusion, lattice-Boltzmann athermal equation, actually, deals with ther- 
mal problems and thermal and athermal lattice-Boltzmann models will be 
here considered using a single approach. 

There are several features that let the lattice Boltzmann, regular-lattice based, 
framework far away from what it would desirable to be its starting point: the 
continuous Boltzmann equation. These features include the particles model, 
the collision model, the polynomial development used for the equilibrium dis- 
tribution, the time discretization and the velocity discretization. Some of these 
main features are discussed in the following. 

Collision model. More widely used lattice-Boltzmann collision models are 
based on a Bhatnagar, Gross and Krook (BGK) relaxation term, [2], firstly 
introduced in the lattice-Boltzmann framework by Qian et al., [3], and Chen 
et al., [4]. Thermal lattice-Boltzmann schemes based on the BGK collision 
model use an increased number of discrete velocities and/or include higher 
order non-linear terms in the equilibrium distribution function ([5], [6], [7]), 
temperature dependent weights, [5], and temperature dependent velocities, 
[12]. BGK single relaxation time collision term restricts the models to con- 
stant Prandtl number. The correct description of fluids and fluid flow requires 
multiple relaxation time models (MRT). A two-parameters model was intro- 
duced by He et a/., [8], using two sets of distributions for the particles number 
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density and for the thermodynamic internal energy, coupled through a viscous 
dissipation term. Full MRT models were firstly introduced in the LBE frame- 
work by d'Humieres, [9], derived from the moments method, by making the 
moments and the velocity spaces isomorphic, [10]. The main criticism to the 
moments method is that the highest order of the kinetic moments possible 
to be correctly described with the LBE equation is limited by the number 
of lattice velocities, [11], and high-order kinetic moments are not correctly 
described when all the b-moments in a b-discrete velocities set are consid- 
ered. In currently produced works dealing with applications of the moments 
method, e.g. [1], the main worry is numerical stability and not the description 
of non-isothermal, multicomponent or immiscible fluids flows, which, effec- 
tively, require additional relaxation parameters with respect to BGK models. 
Dispersion equations are used as constraints for the adjustable parameters 
related to the short wave-length non-hydrodynamic moments and numerical 
stability is assured by buffering these higher frequency moments. This ar- 
tificial shielding can be too dangerous for the complex flow structures that 
naturally appear when the Reynolds number increases. A correct description 
of the growing of these flow structures with the Reynolds number requires to 
increase the lattice dimensionality. 

Time discretization. Most lattice-Boltzmann simulations are based on an 
explicit numerical scheme, with second order time-step, 5, errors. Lattice BGK 
models, have been simulated with implicit numerical schemes, [13], [14], or 
LBE modified explicit numerical schemes, [8], with third order time step, 
0(5 3 ), errors. In spite of the fact that, in athermal models, this truncation 
error can be totally absorbed into the physical viscous term, in thermal mod- 
els errors 0(<5 2 ) seriously affects the viscous heat dissipation term, [8]. 
Velocity discretization. Lattice-Boltzmann method is based on a finite set 
of discrete velocities q and weights 0Ui, judiciously chosen so as to ensure 
isotropy for the even-parity rank velocity tensors and, consequently, for the 
fluid transfer properties. He & Luo,[ll], have directly derived some widely 
used lattices (D2Q9, D2Q6, D2Q7, D3Q27) from the continuous Boltzmann 
equation by discretization of the velocity space, using the Gauss-Hermite and 
Gauss -Radau quadrature of the Maxwellian w-polynomial expansions. Unhap- 
pily, excluding the above mentioned lattices, the discrete velocity sets obtained 
by quadrature do not generate regular lattices. In this sense, Pavlo et al, [12], 
proposed a temperature dependent velocity model based on an octagonal lat- 
tice which is not space-filling but ensures the isotropy of 6 th - rank velocity 
tensors. 

In this work, we present an attempt for deriving the lattice-Boltzmann equa- 
tion, from the continuous Boltzmann equation trying to combine the following 
main features: multiple relaxation-times, 0(5 3 ) time step errors and thermo- 
dynamic consistency in non- isothermal flow. 

In contrast with the moments method, in present conception higher Reynolds 
flows (and non-isothermal flows) require to increase the number of the lattice 
discrete velocities and to increase the accuracy of the LBE equation with re- 
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spect to its continuous counterpart. 

The collision term Q in the linearized Boltzmann equation is modeled by ex- 
panding the distribution function / in Hermite polynomial tensors ^0, which 
forms an orthogonal basis in the Hilbert space H generated by h : C D — > 1Z, 
D being the dimension of the velocity space. Considering that each term C 
(^0 ) is, itself, an element of H, this term is expanded as a linear combina- 
tion of the same order-# Hermite tensors through 2#-order relaxation tensors. 
Isotropy properties are used to reduce these tensors. The infinite series f eq C 
(0) is not truncated. Instead, after a chosen tensor order N, a Gross- Jackson 
procedure is used, [15], and the relaxation tensors are diagonalized. 
It is shown that dHumieres moment equations [10] are particular discrete 
forms of the derived model when the diagonalization constant is considered to 
be zero. 

By performing a Chapman-Enskog analysis of the derived continuous model it 
is shown that the thermo-hydrodynamic macroscopic equations are correctly 
retrieved with a second-order model. Third-order models are only required for 
describing third-order coupling (Soret and Dufour effects) in multi-component 
systems ([16]). 

The derived kinetic model of the continuous Boltzmann equation is then dis- 
cretized. 

It is shown that an explicit numerical scheme with 0(5 3 ) time step errors can 
be derived, using He at al. procedure, [8]. 

Velocity discretization is the most critical step in deriving lattice-Boltzmann 
equations. 

For each N, the equilibrium distribution is taken as an n'^-degree Hermite 
development of the Maxwell-Boltzmann (MB) equilibrium distribution, in ac- 
cordance with the constraints that are imposed by the physical problem. 
Although, as it was shown in [11], velocity discretization of the most widely 
known lattices can be achieved by Gauss-Hermite and related quadratures, 
quadrature schemes did not succeed in generating multi-speed regular lat- 
tices, suitable for thermal problems, placing an, still, open question in the 
LBM framework. In present paper, it is shown that the integrating variable 
has an important role in defining the equilibrium distribution and the lattice- 
Boltzmann model: a) discretization based on the particles velocity c, giving a 
set of discrete, constant, particle velocities cl, leads to temperature dependent 
weights uji (TDW), b) discretization based on the, temperature dependent, 

— * 

dimensionless velocity C gives a set of temperature dependent particle veloc- 
ities cl (TDV). 

In this context, it is shown that the thermal part g e9 in He et al. two-distributions 
model, [8], can be formally retrieved from TDW models, as an w-polynomial 
approximation with errors O (uQ) where 6 is the temperature deviation. 
Although a more complete theoretical analysis is still necessary, the consider- 
ation of TDV models appears to be suitable for thermal problems. However, 
in this case, the particles allocation rules, required for the local conservation 
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of mass, momentum and energy, make the particles number-density, n, to be 
temperature dependent and the implicit temperature dependence of n is dif- 
ficult to manage, when performing a Chapman-Enskog analysis of the kinetic 
models. Finding the macroscopic behavior of these models is, still, in progress. 
A simulation scheme suitable for non-isothermal problems is presented for the 
TDV model. 



1 Boltzmann equation in the continuous velocity space 



1.1 Development of the distribution function in Hermite polynomials 



!• 



- c2 *o, { r 9 )*e, {Se) dC = A e A ( , )( ), (2) 



The Maxwell-Boltzmann equilibrium distribution, [15], can be written as an 
infinite series of Hermite polynomial tensors ^e,(r ), [17], 

where (r$) is a sequence of indexes n, r 2 , ...re and repeated index means sum- 
mation, = 1,*1,q = 2C a , ^ 2 ,Q/3 = 2(C Q C /3 - \5 a p), ^3,a/3 7 = §(C a C/3C 7 - 

^SapC-y — \$a^Cf3 — ^b~f3-yCa) an d so on - The dimensionless particle velocity is 

/ sl/2 

C — ( ^ J c. These tensors are orthogonal in the Hilbert space H, satisfying 

where ^ ^ is a 2#-order isotropic tensor, [18], and \ e is a constant. With 

U = (^f) ^ m, the coefficients a e e q ^ r ^ in Eq. (1) are the moments a e o q = n, 
a eq a = nU ai a e 2 q a/3 = nUJdp, a^ 9 = nUJJpU^ and so on, which are dependents 
on the volumetric number of particles n, on the dimensionless macroscopic 
velocity U and on the temperature T. 

For each point x the distribution function <fi in the non-equilibrium part f neq 
= f eq (j) can be developed in terms of the orthogonal basis ^e,(r ), [16], [19], 

written in terms of the velocity fluctuation C/ = 2 ^~" 1/2 = ^ 1/2 = C — U 

V m ) \ m ) 

^ = E a },(r,)(^*)^,(r,)(^/), (3) 



and coefficients Oq can be related to the macroscopic moments of /. In this 
way, clq = 0, af a = 0. The coefficient a\ a p is related to the viscous stress 
tensor r Q/ g through 

«2,a/3 = (4) 
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where P = nkT is the thermodynamic pressure. 
The fluctuation kinetic energy E(x, t) is given by 



E(x,t) — J f^ 171 (c — u) 2 dc = J f eq -m(c — u) 2 dc. 



In this way 



or 



In two-dimensions 



or 



f neq \m (C) 2 dC = 0, 



I f neql -mC a C a dC =Ur (r) = 0. 



a 2,xx + a 2,yy ~ 0- 



For third-order moments 



S a /3^ = J fmc a c /3 c y dc = J f eq mc a c /3 c y dc + J f^mCaC^c^dc 



qeq , qneq 

o aj3l -+- o a/3 ^, 



with 



For the non-equilibrium part, 



(5) 

(6) 

(7) 

(8) 
(9) 



(10) 
(11) 



^7 = / f neq mc fa c fl3 c h dc+ (r^u a + T ai u p + r af3 u y ) , (12) 

resulting, using af a = 0, the invariance property with respect to index per- 
mutation and Eq. (11): 



p( 2kT Yn* - Sa * 



m 



2 



\pu oi u0a 1 + \P (5/3 y u a + 5 ai Ufs + 5 af3 u^) 



(13) 



When (3 and 7 are contracted, defining e a to be the total energy flux along 
the direction a, 



I 2kT\ 1 + 



1 



D 



-f)U U a + P y— + l)U a + TapUp 



Qa, (14) 



where q a is the net heat flux along the direction a, i.e., the total energy flux 
e a , subtracting from it, the flow of macroscopic kinetic energy ^pu 2 u a , the 
compression work P + u a and the viscous work r a 



/3U/3. 
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1.2 Collision term 



Particles are supposed to be material points without volume and only able 
to exchange translational kinetic energy, but the collision term Q is, here, 
considered to take multiparticles collisions into account. Since, in this case, 
the collision term structure is not known, some assumptions are required. In 
this manner, near the equilibrium, Q is considered to be f eg £(<f)), the operator 
C being a linear operator. This property was shown to be true for binary 
collisions, [14] and is, here, extended for multiparticles collisions. When / is 
near f eq , Boltzmann equation reads 

d t f + c.vf = n = rc( ( j>). (15) 

Using the development, Eq. (3) , 

^)=E<(r^K(r,))- (16) 

e 

The #-order tensor C (^e,(r g )) is, itself, an element of the C D space and can 
be developed in terms of the 6>-order Hermite tensors that belong to the 
orthogonal basis of this space, 

where 7( rfl ),(s fl ) designate the (rg), (sg) components of 2#-order relaxation ten- 
sors. Considering, as for binary collision, £ to be a self-adjoint operator, with 
non-positive eigenvalues, 

je~ c2 fC(^ B( .W fl , ,dC f 

Using Einsteins notation 

W) = Y,%.U.Xlr.)*..<..V ( 19 ) 
6 

where repeated indexes mean summation. 

Above equation is an infinite summation on 0. When the terms above a chosen 
order N are diagonalised, following a Gross- Jackson procedure, [15], 

N oo 

£ W = EW<M^ (S9l -7 N+I E S (r e) , (se) a i(r e )*e, (se y ( 20 ) 
0=o e=N+i 

where 

S, w . — 5 r 8 ....5 r 8 . (21) 
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In this way, using Eq. (3) 



N 



£(</>) = 



u>=o 



(22) 



where A, w , = — (7, w , +7„, 1 5, w J is positive for all ro,sa, since a) 
A (r } (s } = — 7 ( ) (s j for all off-diagonal components and b) the diagonal com- 
ponents 7 (r , ) (r ) are negative with an absolute value that is greater than 7 Ar+1 
for all 9 smaller or equal to N. Eq. (22) can be considered as an N^-order ki- 
netic model to the collision term, with an absorption term / y N (p resulting from 
the diagonalization of the relaxation tensors after the given N. Therefore, all 
the moments of order higher than N are collapsed into a single non-equilibrium 
term minimizing the truncation effects on the fine structure of the operator C 
spectrum. 

Although very little is known about the true collision term fl when multiple 
collisions are considered, Eq. (22) generates increasing accuracy models to Q 
when the distribution function / is near the Maxwell-Boltzmann equilibrium 
distribution, f eq . The only restrictions are: a) particles were considered as ma- 
terial points without volume and b) particles internal energy and long-range 
forces among the particles were not considered in the derivation. 
When N — or N — 1, Eq. (22) gives the well known BGK model, when all 
the collision operator spectra is replaced by a single relaxation term. 
Each term in the sum, in Eq. (22), gives the relaxation to the equilibrium of 
second or higher order kinetic moments Mg that are not preserved in collisions, 
modulated by a \g relaxation tensor. In Section 1.3, explicit expressions are 
given for the collision models. When the diagonalization constant is consid- 
ered to be zero, i.e., when the series, Eq. (19), is truncated above N, replacing 
q _ y e ?£(0) j n the Boltzmann equation, the inner products of the resulting 
equation by s S> x ( } give 

d * a L(r x ) + ^ a L(r x ) = \r x )(s x) a n x l x) , (23) 

where the distribution function / was developed following 

/ = E e " C? <(r fl) *e,(r 9 )> ( 24 ) 
e 

and ^ 

a ^) =n (^) <<'•>■ (25) 

It can be easily seen that DHumieres moment equations ([9], [10]) are partic- 
ular discrete forms of Eq. (23). Nevertheless, in dHumieres moments method, 
all the b-moments in a b-discrete velocities set are considered. It was shown 
in [11], that the number of degrees of freedom of a given lattice restricts the 
order n of the kinetic moments with exact quadrature. This means that all 
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the moments which order are greater than n cannot be correctly described 
in this given lattice. As it was mentioned in the Introduction, in the mo- 
ments method these high-frequency moments are forced to give consistent and 
numerical stable low-frequency macroscopic equations by using dispersion re- 
lations, decreasing the effect of numerical instability sources, but buffering the 
appearance of complex flow structures, when the Reynolds number increases. 



1.3 Collision models for the continuous Boltzmann equation 



In present section, the isotropy of A th and 6 th rank tensors will be used to 
give explicit forms for the second and third-order collision models, Eq. (22). 
Without any loss in the generality, we restrict ourselves to two-dimensional 
spaces. 



1.3.1 Second order model in the two-dimensional space 
From Eq. (22) 

\r 2 ),(. 2 ) a 2,(r2)*a,(. a ) = KfhS^afi^W 



(26) 



Requiring isotropy of 4 th rank tensors and considering the symmetry with 
respect to index permutation, 



(27) 



In this way, 



K/3-yS a 2,a/3^2 n 6 - -V °2,aa^2, 77 + «2,a/3^ ' 2, a f3 + a 2,a/3* ' 2,0 a 



A, 



a i,xx {f-}x 2) + a tyy { C fy 2 

2 a 2, X yCfxCfy 



(28) 



since a\ a(x = 0. Using Eq.(4) 



T~xx I Cf x 



H~ T~vv I C 



yy \ ^fy 



2 T xyCfxCfy 



, (29) 
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or, from Eq. (8), r xx = -r. 



mi 



Txx {f- fx Cfy) + 2T xy Cf x Cfy 



(30) 



Second order model in two dimensions will be written as 



£ (2) (0) = "^ 



^TxyCfxCfy 



(31) 



Third- order model 



From Eq. (22) 



(32) 



For isotropic fluids, tensor A a/374C)7 is a linear combination of five 6 order 
tensors given by the recurrence relation, [18], 



A (6) = g A (4) 



(33) 



resulting 



0/3 



(34) 



1.4 Macroscopic thermohydrodynamic equations 



Macroscopic thermohydrodynamic equations may be obtained from the Boltz- 
mann equation by multiplying this equation by the mass m, the momentum 
mc and the kinetic energy |mc 2 of the particles and integrating the resulting 
equations in the c velocity space. 
The mass conservation reads, as usually, 



d t p + d a (pu a ) = 0. 



(35) 
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/,From the momentum preservation in collisions 



d t {pu a ) + d a {pu a U(3 + P5 a p + r a p) = 0, (36) 

where P is the thermodynamics pressure, P = nkT and r a p is the viscous 
stress tensor. 

When Eq.(36) is multiplied by u a , the macroscopic kinetic energy, \pu 2 ) bal- 
ance equation is obtained 



d t ( K \p U l) = PV - U + TapdpUoc - 9/3 Q^pulup + P5 a pU a + T a pU^j . (37) 



The total energy conservation equation reads 



dt (e + ^P«a) = -fy {^P u l + E^ju (3 + (P5 a p + T aP ) u a + qp 
where E is the thermodynamics internal energy 



(38) 



E= l\m (c - uf fdc= f \m (c - uf f eq dc= ^nkT. (39) 



The internal energy balance equation is obtained by subtracting Eq. (37) from 
Eq. (38), 



d t (E) = - (Py.u + r a pdpu a ) - dp [Eup + qp] , (40) 
where — (PV.-u + r a pdpu a ) is the source term of internal energy. 

Equations (35, 36 and 40) form a closed set of equations when the viscous 
stress tensor r a p and the heat flux vector qp are known in terms of the spa- 
tial gradients of the first macroscopic moments, p, u and T of the distribution 
function. This is accomplished when the Knudsen number, Kn — > 0, by per- 
forming a Chapman-Enskog asymptotic analysis of the modelled Boltzmann 
equation. 
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1.5 Chapman Enskog analysis for the continuous model 



Considering f° in the asymptotic expansion 

f = f + Knf 1 + (41) 

to be the Maxwell-Boltzmann equilibrium distribution f eq (n, u, T), the zeroth 
order time derivative resulting from the Chapman-Enskog induced decompo- 
sition of the time derivative reads, 



1 d f° 
f° dt 



+ 



(2kT\ 



1/2 



y m 



(42) 



1.5.1 Second order model in two dimensions 



Using Eqs. (42 and 31) 



2 (c% - ^ d x u x + 2 (c% - ^ d y u y + 2C fx C fy {d x u y + d y u x ) 

V-+(^) 1/2 (^-2)C>lnT 



C % 2 



C /y- 2 



2A, 



7 1 C% 9 ) + T yy [ C% 2 ) ~^ ^ Tx y^f x ^fy 



(43) 



For finding the correct expression of T a @, in terms of the spatial derivatives 
of the macroscopic variables, the inner product of the above equation by 
(Cf a Cfp — \b~ a p) in the C/ velocity space is performed. By multiplying the 

above equation by (Cj x — |) 



T~xx 1 1 

( A M + 7s) = -drlla + - V.U = - (d y Uy - d x U x ) 



(44) 



Similarly 



T 1 

( a m + 7 3 ) jr = 2 ( 9 * u * - 3/%) ' 



(45) 
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and 



(^ + %) 1 f = -\(9 x u y + a y u x ). (46) 
These results give for the first and second viscosity coefficients, 

2 \ + 7 3 

in the relation 

r a /3 = -fi (d a u p + dpu a ) + rjs/.u. (48) 

Eq. (47) means that the first and the second viscosity coefficients are not 
independent quantities and this result must be considered as a limitation 
resulting from the continuous collision term itself, where the particles were 
considered as material points with translational degrees of freedom. In fact, 
this same result will be retrieved when using the third or higher order model 
for the collision term. The consideration of internal energy modes would be 
necessary for an up-grade of Eq. (47), [16]. 

The third-order moment at xxx may be obtained by multiplying Eq. (43) by 
(C% - f)C fx . In this way, 



3 \ 777. 



1/2 



at, xxx = -j^-^d x \nT. (49) 



Multiplying Eq. (43) by (C% - ±) C fx , 



1 



/ 2kT\ 



1/2 



<yy X = -J^Z^dJnT. (50) 

^ /3 



But 



«* = P (^) V2 (<C + <vv*) = ~^T d * T > ( 51 ) 
giving for the thermal conductivity 



{D + 2 )n ^T L 

2m % 1 ; 
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In this manner, present second-order continuous kinetic model is thermody- 
namic consistent and able for analyzing non-isothermal and fully compressible 
flows. The thermal conductivity is related to 7 3 diagonalization constant. Con- 
sideration of third-order models will be, only, necessary in multi-component 
systems, for correctly describing third-order coupling: the Soret and Dufour 
effects, [16]. 



2 Discretization 



In present sections an analysis is performed, trying to emphasize the theo- 
retically identifiable effects of time and velocity discretization on the ability 
of the derived discrete models in retrieving the correct thermohydrodynamic 
equations, i.e., the full compressible Navier-Stokes equations and the thermo- 
dynamic internal energy balance equation, with the Fourier heat flux term. 



2.1 Time discretization 



Boltzmann equation with the kinetic model Eq. (22)becomes 



d N 

jj + ij = r E K e usAir-^^ e) + iJ eq - (53) 

For avoiding time-step errors 0(5 2 ), Boltzmann equation is integrated be- 
tween t and t + S, considering linear approximations for f eq (x + ct', c,t + t') 
and, also, for a^ r } , since a^ r ) (x + ct' ,t + 1'), when < t' < 5, [8]. The 
result is 



/ (x + c5, c, t + 5) - f (x, c, t) 



: (InS) - [r + 85, c, t + 8) + r (x, c, t)} 
- (InS) \[f(x + cS, c,t + 6) + f (x, c, t)} 

N 

+ EK,, (s > 9 ,„,/ e? (v^) x 

X 2 K(^) ^ +cS,t + 5) + aj >(rfl) (x, t) 



(54) 



This corresponds to an implicit numerical scheme (in fact, a Crank-Nicholson 
scheme). Although implicit schemes are very easily manageable in the lattice- 
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Boltzmann context, [13], [14], if one wants to avoid implicitness a new distri- 
bution can be defined as, 



1 N fi 

f = f + ^Hf- f eq ) - E 2 \r e UsX(r e )*e, iSe) f eq , (55) 

8=0 



resulting 



/ (x + c5, c,t + 5) 

'at i n U 



N 1 2 

s N 

-—^1$ S T "\r e Us e ) a t(r e ) f^e,^ ) > (56) 



where = It must be observed that at (r ) are the macroscopic moments 
of / and not of /. Nevertheless, it can be shown from Eq. (55) that g$( s ) and 
a t (r ) are directly related by 



Consider, for instance, the second order model. In this case, it can be shown 
that 



= -2A„ (a^) , (58) 
since a^aa = an d V is required to be positive. Eq. (57) means 



~at rf6 =[l + — + 5\^al fS . (59) 



2.2 Velocity discretization 



Velocity discretization is the most critical step in presently proposed proce- 
dure. For athermal problems He and Luo, [11], have shown that some widely 
used sets of discrete velocities = 1, ...&} may be derived from the con- 
tinuous velocity space by Gauss-Hermite (D2Q9, D3Q27) and Gauss-Radau 
(D2Q7) quadrature. All these sets are space-filling, in the sense that for every 
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lattice site x , x + c\ points to another site in the lattice. 
If it is agreed that quadrature is the bridge connecting the continuous and the 
discrete velocity space, discretization means to replace the entire continuous 
velocity space c D by some discrete velocities satisfying the quadrature for 
all the kinetic moments of interest, i.e., for all the kinetic moments that are to 
be correctly described in lattice-Boltzmann simulation. Although it is highly 
desirable set q to be space-filling, this condition is not essential for the dis- 
cretization itself. 

When performing the quadrature, an integration variable must be chosen. If 
the dimensionless fluctuation velocity Cf = - 2 ^y /2 is chosen as the integrat- 
ing variable, considering % to be a polynomial of degree r in the velocity, 



< 



X >= f f X dc = n^— f e- C f X ' (<?/) dC f = n £ cu iX ' (Cfi) , (60) 

where Cfi is a discrete velocity (a constant vector), dependent, basically, on 
b and on the kind of quadrature operation it is being performed, X ' (pf) * s 

a polynomial in Cf of degree r p = r +s when it is related to a preserved 
moment and r p = r+s+1, otherwise, [11], s being the degree of the polynomial 
approximation to f eq and Ui are the constant weights to be attributed to each 

— * 

discrete velocity Cfi. Exact quadrature restricts the highest value of r p to r Pm . 
For a given class of quadrature we can write r Pm = r Pm (b), in the sense that 
increasing b enables higher degree polynomials to have exact quadrature. 
For the first kinetic moment, n, 



n =< 1 >= n— ^ [ e c hdC f = n^u^l = n^cu,, (61) 
71 J i=i i=i 

resulting, 



if = coin. (62) 

This means that the discrete equilibrium distribution does not depend, explic- 
itly, on the macroscopic velocity u and on the temperature T. These depen- 
dences are included in the particle velocities through, 



(2kT\ 1/2 - 

Ci = u+l— -J C fi = Ci{T,u). (63) 

When the Cfi generate a regular lattice, this choice is possible in LB frame- 
work, since LB simulation may, at least in principle, be performed in the 
(x,Cfi) space. Nevertheless, the physical grid (x, q), i.e., the physical grid 
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points where the particles will be located after each time step, will be time de- 
pendent, simulation tends to be very cumbersome and, at a first sight, bound- 
ary conditions will be difficult to satisfy. 

Another choice is the dimensionless particle velocity C = ^ 2k £y/2 ■ This re- 
quires to rewrite the equilibrium distribution as in Eq. (1), but, now, the 
series 

E a ?(r,)M>m*,(r.)> (64) 

e 

must be truncated somewhere. This is an important distinguishing point of 
discrete models, since in each continuous model presented in Section 1.3, al- 
though the full collision term is replaced by its N i/l -order approximation, the 
equilibrium distribution is, always, the full MB distribution. 
Second-order approximations are widely used in athermal simulation, but ther- 
mohydrodynamics require third (or higher) order approximations for the equi- 
librium distribution, 



1 + 2C a U a + 2 [CaC/3 — \5 a p) U a U(i+ 

4 / 3 \ ' ( ' 

3 [C a C/3 — 2&af3) C^UJAplA^ 

which can be viewed as a third-degree polynomial expansion of the f eq depen- 
dence on U, with errors 0(W 4 ). 

After quadrature, the equilibrium distribution becomes 



f eq = n 



ir D / 2 \2kT 



m 



„ ea I 1 + 2CiJJ a + 2{CiaCip — \5 a p)U a Up 

— * 

where, as above, the weights Ui and the velocity vectors Ci are dependent on 
b and on the kind of quadrature that was performed. 

When u — 0, the equilibrium distribution is only dependent on the tempera- 
ture T through the number density of particles, n. Nevertheless, the particle 
velocities are temperature dependent, 




c1=^J C t = cdT). (67) 

A simulation alternative is presented, in this case, by redistributing the par- 
ticles among adjacent sites, in accordance with allocation rules, locally pre- 
serving the mass, momentum and kinetic energy of the original packet. This 
strategy will be discussed in Section 2.3, leading to the establishment of tem- 
perature dependent velocity models (TDV). 
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Pavlo et ai, [12], developed a TDV model based on an octagonal lattice, 
which is not space-filling but assures the isotropy of 6 th rank tensors. It can 
be shown that this octagonal discrete velocities set can be retrieved using a 
Gauss-Radau quadrature, with 8 angular directions, giving r Pm = 7, instead 
of 5 as in the D2Q7 model and assuring the exact quadrature of third-order 
moments. Additional considerations can be found in Section 2.4. 
Avoiding the q temperature dependence requires to consider the particles ve- 
locity c as the integrating variable when performing the quadrature, i.e., to 
let c 2 free from T in the exponential part e _c of the equilibrium distribution. 
This can be accomplished by writing 



e^F = ( e - c !o) - > ( 68 ) 
where T is a reference (and constant) temperature and Cf Q = 2fc c T ~" 1/2 is a 

( m ) 

new dimensionless fluctuation velocity referred to the temperature T Q . When 
T is near T , i.e., when the departures from thermal equilibrium are small, 
the above expression may be developed in a Taylor series around — 1. 
Considering Q — -S- — 1, this development gives 



-c? 
e f° 



i + c%e+ 1 -c%(c%-2) e 2 + ... 



(69) 



which terms are increasing powers of Cj . 

In this way, retaining just the first power in 9, 



D/2 



1 + 2Co iQ Wo, a + 2 (Co, a Co,/3 ~ Mq,cMo,i3 
+ | (CoaCo/3 — §<W) CoryUoaMopMo'y 



x e 



(70) 



where Uq 



I 2kT Q \ ' 
V m / 



In this case, the quadrature will give for the discrete equilibrium distribution, 



fT = 9i[T,C 



1 + 2Co,i a Mo,a + 2 (Co,i a Co,if3 ~ §<W) ^0,a^0,/3 



, (71) 
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where 



,rp s D/2 

9t (T,C%,,) = (^) [l+C^e], (72) 



is a temperature dependent weight. When T > To, gi reduces the amount 
of particles with zero velocity, redistributing them to the kinetic modes % in 
accordance with C| 0i and inversely when T < T . It is well known by LB 
practitioners that this redistribution is highly desirable in LB simulation and 
redistribution rules were empirically found by some authors (e.g., [5]). 
Considering a linear approximation to the temperature non-equilibrium, 



T \ D ' 2 D 



T 

Eq. (72) can also be written as 



) « 1 " y6. (73) 



9i (t, c% :t ) = i - y e + cj 0ii e + o (e 2 ) , (74) 

and the equilibrium distribution may be written as a sum of two distributions 



f? = f?n+ nh (75) 

where, dropping-out the third-order term, 



fi,T = ®C%,i^in 



1 + 2Co t iaUo t0l + 2 (Co t i a Co,i/3 ~ ^afij ^0,J^0,/3 ■ (76) 



This equilibrium distribution is related to the thermal distribution function 
g in He et al. two-distributions model, [8]. To fit the model into a D2Q9 

lattice, He et al. have, further, replaced Cj 0i by [Cq^ — Uoj , truncating all 

the terms O an( i higher, after the multiplication. It can be easily seen 

that the resulting expression has second-order errors O (©W ) limiting He et 
aVs model to low local speeds. 

We have preferred a somewhat different decomposition in Eq. (68), working 
with the particles velocity c and not with the fluctuation velocity (c — u), 
making 



.. 2b 

2kT 

e m 



= (e- C °) T , (77) 
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resulting in a temperature dependent weights model (TDW), which equilib- 
rium distribution is given by 



n 



eq 



9i 



UJiU 



1 -+- Z h Z ^Lo,ia'-'0,i/3 ~ 2 To a /3 J ~ ^~ 3_ 



+ 3 \y0iaL>0i(3 ~ 2 7\) a/3 ) 



(*)' 



(*)" 



, (78) 



where 



^(T,C 2 l )=l + (c 2 i -^) 0. 



(79) 



Since |r, (|r) , ... in Eq. (78), appears inside the polynomial expansion these 
terms must be, also, developed in 6 for preserving consistency in the order of 
approximation. After multiplication, an expression for the equilibrium distri- 
bution in 2D problems can be written as 



fi =Uin< 



C 2 



1)6 + 2 



C 2 



2 



+2(C ,i.Uo) -Ui + l(C ,i.Uo) -2(C ,i.Uo)U 2 



(80) 



with errors O (0 2 , U%Q). A redistribution expression similar to the above one 
was obtained by Shan and He, [19]. Large temperature deviations require to 
consider additional powers in 0, in the Taylor expansion, Eq. (79) and to 
increase the number of discrete velocities in the lattice. This will be discussed 
in Section 2.4.3. 



2.3 TDV model: allocation rules and collision step 



Considering the analysis shown in the last section, TDV models appear to be a 
promissing alternative for LB simulation of non-isothermal problems, since no 
theoretical limitations were identified, related to the thermal non-equilibrium 
deviation 0, as in the TDW models, Eq. (80), where the Taylor series in 0, Eq. 
(79) was truncated after the first-order term. In present section this model is 
discussed, considering that temperature-dependent velocities require to mod- 
ify the propagation step in the LB simulation scheme. 

After quadrature, the third-order approximation to the equilibrium distribu- 
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tion Eq. (66) can be written in a general form as 



1 + 2^^ + 2^(0^-^^ 



7 



+ 3° (, C fca C fc/3 - 4^%°^ 



(81) 



where a is a lattice constant related to the lattice symmetry. The u k are the 
lattice-weights, defining the inner product 



f-9 = ^2^kfk9k- (82) 

k 

In general, the discrete velocities set generated by quadrature will be not space- 
filling requiring simulation strategies similar to the one proposed by Pavlo et 
al, [12]. In space-filling lattices such as, e.g., the D2Q13H, the equilibrium 
distribution will keep the form given by Eq. (81) but the weights u k must be, 
empirically chosen for ensuring isotropy of even-parity rank velocity tensors 
up to the 6 t/l -rank. 
Particle velocities are given by 



c k = - 5 ? k - (83) 



where the lattice dimension h is given by 



h 
1 



'2kT 



m 



1/2 



(84) 



and 



cl = W^4, (85) 
V J o 

where {e k , k — 0, b} are the lattice unity- vectors. 

In present scheme, macroscopic velocity u* = uj (jf), is calculated as nu* = 
J2k /fcCfe where n = J2k fk- Temperature is found to be such that 

n^=a 2 J2f k (cl-u*) 2 . (86) 



After discretization, the second-order model reads, in two-dimensions reads 
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x6a 2 Ta_!—f<" 



Ms) 



2-/fc 



'xx \ C fk,x 2a? T ) ^ 
W \ fk,y 2a 2 T ) ~ xy^fk,x u fk,y _ 



where r 3 = 1/% and r M = 1/A„. 

The dimensionless viscous stress tensor is calculated as, 



Taf3 



'a/3 



2a 2 nkT 



1 V~~"^ rneq * * 
— 2-^Jk c ka c kf3i 
71 k 



and the heat flux, as 



(87) 



q 1 



Erneq * *2 
Jfc C fk 1 C fk- 



(89) 



The modified fluxes r*^ and q* needed in the simulation have similar expres- 
sions in terms of f^ eq . 

The particles that are present at site x after the collision step have a ve- 
locity S* = When < 1, the kinetic energy these particles have 
is just enough to enable them to jump to some intermediate position be- 
tween x and x + Si. In present discrete model these particles are redistributed 
to the next contiguous sites x and x + Si preserving the mass, kinetic en- 
ergy and the momentum of the original particles packet fi (x, t). This is per- 
formed by using a levers rule, making f i0 (x,t + 5) = fi (x, t) (l — and 

fi _i (x + Si, t + 5) — fi (x, t) (ifr) , and by imposing to both amounts /^o and 

fi _i the same non-integer velocity S* related to the particles-packet before 
the propagation. In this way model particles are velocity-memory particles in 
the propagation step. In the notation fy, the index % is related to the site 
direction and j means from what site these particles were come: j = means 
that the particles were come from the same site, j = —1, from the site x — Si, 
j = —2, from the site x — 2e*j and so on. 

In this manner, the particles propagated along the direction % are redistributed 
among the departure site and the next-one contiguous site, but have their ve- 
locity unaltered, after the propagation. These rules preserve, locally, mass, 
momentum and energy and non-physical convection is avoided. 
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When 1 < < 2 , particles will be redistributed following 



/;,_! (x + e u t + 5) = ft (x,t) i ^4 - ) 



(90) 



and 



/<,-2 (f+2e i ,* + <J) = / i (f,t) 



1 / T 



3 VTr 



- 1 



(91) 



Temperature ratios y^T greater than 2 are not expected in present work, 
since they would represent very strong thermal non-equilibrium. 
After propagation, at time step t, all the distributions f it j (x, t) are known in 
each site x. The amounts fij (x, t) enable to calculate the equilibrium moments 
n(x,t), u*(x,t) and T(x,t) and the equilibrium distribution, Eq. (81), with 
k = (i,j). 

The non-equilibrium distributions are given by f™J q = fij — f if and the viscous 
stress tensor can be thus calculated. In the collision step, collision will give a 
single distribution fi (x, t) from the several j amounts fij in the direction i 
of the site x at time t, 



hi + (fij fij) 



6a 2 H 



2t 



N 



S 

' 2t u 



1 feq 



( c f. . _ 1 

V f,*3,x 2a?T J 1 'yy 

-4-2t* r* r* 

~ xy f,ij,x f,ij,y 



X 



T \ 1 ~* ( *2 _ 1 T_ 
r /„„ y-f,ij, y 2a 2 " T 



(92) 



These particles will move to the next sites with the site velocity c*i — , , 
where T = T(x, t) is the site temperature, calculated using Eq.( 86), just after 
the propagation step and with k replaced by 

In the propagation step, a fraction of the low-speed particle-packets, in sites 
where T is small, will remain in its departure site. The overall effect is to 
increase the number-density of particles in the sites where the temperature 
is small and to decrease it in the sites where the temperature is high, emu- 
lating the temperature dependence of the MB distribution. In this way, the 
equilibrium distribution is temperature dependent even when the macroscopic 
velocity u is zero, i.e., when f^ 9 = UiU, because n is temperature dependent. 
Also, the temperature derivative of the equilibrium distribution will vary in- 
versely with the temperature, similarly to its continuous counterpart. This is 
an important condition for retrieving the temperature derivative term in the 
second member of Eq.(42). 

Nevertheless, the temperature dependence of n is related to the particles al- 
location rules used and further theoretical analysis of presently proposed pro- 
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cedure leading to a correct Chapman-Enskog analysis of TDV model is, cur- 
rently, under investigation. 



2.4 The nine bits lattice 

In this section, the TDW and TDV models are explicitly written for the D2Q9 
lattice and their feasibility for simulating isothermal and non-isothermal prob- 
lems in this lattice is discussed. 



2.4.1 TDW model 



The D2Q9 lattice is not suitable for thermal problems since: i) 6 th rank ten- 
sors are not isotropic, ii) energy transfer q a ^ is not correctly described in 
this lattice and iii) the lattice-dimensionality is too small. However, the use of 
thermal models in the D2Q9 lattice can, possibly, be shown to be a promising 
alternative for simulating near-isothermal problems, when the main interest 
is to increase numerical stability, with respect to athermal simulation. In fact, 
in its kinetic theory concept, temperature is varying from site to site, with the 
local velocity u. When this variation is not considered, i.e., when the growing 
of temperature gradients are not annihilated by heat flow, they can, possi- 
bly, become instability sources in the numerical scheme. Actually, in contrast 
with classical CFD simulation of isothermal Navier-Stokes equations, lattice- 
Boltzmann simulation is never isothermal and the consideration of heat flow 
should be helpful in the simulations. 

For the D2Q9 lattice, r Pm {b) = 5 = 2X3 — 1. Gauss-Hermite quadrature gives, 



where e t = (0, 0), (1, 0), (0, 1), (-1, 0), (0, -1), (1, 1), (-1, 1), (-1, -1), (1, -1). 



The weights Ui are found to be u = 4/9, Ui = 1/9, i = 1, ...,4 for the main 
axes and Ui = 1/36, i = 5, ...,8, for the diagonals. 

Using a dimensionless macroscopic velocity u* = x and dropping-out the 



third-order term, since, as commented above, the exact description of energy 
transfer is out of thought in this lattice, the equilibrium distribution in the 



[11] 




(93) 



5 
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TDW model will read 



= UsTl 



i + l( e f-|)e + 3[i + l(ef-|)e 



+ | (ei.it) 



2 
3 



U ' 



(94) 



It can be easily seen that the equilibrium distribution given by Eq. (94) satisfies 



n 



(95) 



(96) 



P = nkT. 



(97) 



The internal energy 



the momentum flux 



E = J2f! q \m(c l -ur, 



(98) 



(99) 



and the heat flux 



q e a q ^Efi 9 om(&i-u) 2 (c ia -u a ) : 



(100) 



are also retrieved as, respectively, ^nkT, P5 a p + pu a up and 0, with errors 

o(w 2 e,e 2 ). 



2.4.2 TDV model 

All the restrictions given by Eqs. (95-97 and 98-100) are satisfied when the 
equilibrium distribution is given by Eq. (81), dropping-out the third-order 
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term, i.e., 



f? = u> k n l + 3^ + l(cl a cl $ -~S aP 



) 




(101) 



Particles redistribution in TDV model should increase numerical stability with 
respect to the athermal simulation method, since, in this last scheme, particle 
allocation is not temperature dependent and sites with higher T will rest with 
an excess amount in the number of particles, contributing to the enhancement 
of non-physical source-terms. 

In Pavlo et al. , [12], TDV model, the use of an octagonal lattice, instead 
of the D2Q9 lattice, was considered for ensuring isotropy for the 6^-rank 
tensors. A second-order interpolation scheme is used for allocating the particles 
among adjacent sites in the propagation step, following a procedure that is 
similar to the one exposed in Section 2.3, but the model particles are not 
velocity-memory particles and the D2Q9 e*j lattice- velocities are attributed to 
the redistributed packets. This is performed preserving mass, momentum and 
energy. After the propagation, the number density, the momentum and the 
kinetic energy are calculated in each site. Previously to the collision step, the 
distributions /, are renormalized to account for the octagonal lattice velocities, 
satisfying the new local moment constraints. Successive transitions between 
the D2Q9 and the octagonal lattice can affect the isotropy properties of the 
6^-rank tensors and, in our opinion, this is the only question that remains 
to be answered, since numerical viscosity effects can be avoided by reducing 
the spatial scale. In this manner, the problem is, nearly, the same as above, 
because the correct description of the heat flow vector cannot be assured in 
this model. 



2.4-3 Increasing the number of discrete velocities 

After the very promising results of He and Luo's work, [11], when the D2Q9, 
D2Q7 and D3Q27 lattices were found by exact quadrature, it is, now, appar- 
ent that the integer multiplicative factors that are required by regular-lattice 
velocity components are a too strong restriction: for their exact quadrature, 
Maxwellian-shaped curves require the roots C{ a to be progressively close when 
their absolute values increase. 

For some regular lattices, the weights uJi can be found to give even-parity 
isotropic tensors. In this manner, isotropy up to the 6* A -rank tensors is as- 
sured for the D2Q13H lattice by choosing ujq = 132/300, for the zero velocity, 
= 27/300 for the first 6 velocity vectors, with |e*j| = 1 and 002 = 1/300 for 



the next 6 velocity vectors, with |e*j| = In this lattice h/S — y 10 ^ ■ 

Simulation of non-isothermal problems is possible in this lattice by using third- 
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order approximations to the MB equilibrium distribution, derived, alterna- 
tively, from Eq. (80) in the TDW model or from Eq. (81) in the TDV model. 
In TDW models, improving the approximation, Eq. (80), to the MB equilib- 
rium distribution, considering 9 2 powers in Eq. (69), also requires to take 
4^-degree powers of C into account and, consequently, to consider 4 th - de- 
gree polynomial expansions in Uq for the equilibrium distribution. This has 
been performed for a D2Q21S lattice, choosing the weights u>i for giving even- 
parity isotropic tensors up to the 6* h -rank and for retrieving the equilibrium 
conditions, Eqs.(95-97 and 98-100), with errors O (W 3 6, 9 3 ,W 5 ,W © 2 , •••)• Th e 
D2Q21S is composed by the zero velocity and by the superposition of five 
square lattices with a 45° shift between each two sequential lattices, with 
speeds 0, 1, y/2, 2, 2\[2 and 3. The equilibrium distribution was found as 

1 + (ef - 1) + (1 - 2ef + ef/2) 9 2 + 2 [1 + (e 2 - 2) 9] e t .u*~ 
+2 (e-.u*) 2 (1 + (e 2 - 3) 6) - (u*) 2 (1 + (e 2 - 2) 0) 
+ | {ei.u* f - 2 (ei.u*) (u*) 2 + § (e^*) 4 
-2(e i .u*) 2 (u*) 2 + \(u*f 

(102) 

with TOpio-ko / i . — / J79_ _41_ _5_ _31 1 1_\ 

wiui weignib w t — | 1152 , 384 , 96 , 3840 , 1536 , 5760 J- 

3 Conclusion 

As a kinetic method, LBM is based on a discrete and finite approximation to 
the continuous Boltzmann equation giving an approximated solution to the 
particles distribution function. In this manner, all the macroscopic moments 
of interest will be affected by the accuracy of the modelled LBE when com- 
pared to its full continuous counterpart. 

MRT increasing accuracy models to the collision term f2 were derived col- 
lapsing the higher-order terms, related to the relaxation of the higher-order 
moments, into a single non-equilibrium term, minimizing the truncation ef- 
fects on the fine structure of the collision operator spectrum. 
Thus, in contrast with the moments method, only the hydrodynamic moments 
that are required to be correctly described, are considered in the present mod- 
els. This strategy avoids the use of dispersion relations for buffering undesir- 
able effects of the high-frequency, non-hydrodynamic moments, but claims to 
increase the lattice dimensionality. 

Time discretization requires implicit or modified explicit numerical schemes 
for avoiding 0{5 2 ) time step errors and non-physical source terms in the in- 



f- q = UJiU 
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ternal energy balance equation. 

For the discretization of the velocity space, the requirements for the Boltzmann 
equation quadrature to be exact establishes the minimal discrete velocity set 
that is needed for the lattice, in accordance with the order of approximation 
for the kinetic model that is intended to be used. Nevertheless, these sets do 
not generate regular lattices when multi-speed models appropriated for ther- 
mal problems are considered. 

When performing the quadrature, it was shown that the integrating variable 
has an important role in defining the equilibrium distribution and the lattice- 
Boltzmann model, leading, alternatively, to temperature dependent velocities 
(TDV) and to temperature dependent weights (TDW) lattice-Boltzmann mod- 
els. 

In TDV models no theoretical limitations were identified, related to the ther- 
mal non-equilibrium deviation 6, as in the TDW models. 
Finding the macroscopic behavior of TDV models through a rigorous Chapman- 
Enskog analysis is difficult due to the implicit number density of particles 
dependence on the temperature and is, still, in progress. 
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